module nonlin_solve use iso_fortran_env use nonlin_error_handling use nonlin_multi_eqn_mult_var use nonlin_single_var use nonlin_linesearch use nonlin_helper use nonlin_types use linalg, only : qr_factor, form_qr, qr_rank1_update, lu_factor, & rank1_update, mtx_mult, recip_mult_array, solve_triangular_system, & solve_lu implicit none private public :: line_search_solver public :: quasi_newton_solver public :: newton_solver public :: brent_solver public :: newton_1var_solver type, abstract, extends(equation_solver) :: line_search_solver !! A base class for nonlinear solvers that improve convergence by !! combining a search direction with a line search. At each iteration !! the step is taken as !! $$ x_{k+1} = x_k + \alpha_k p_k $$ !! where the step length $\alpha_k$ is chosen to satisfy a sufficient !! decrease condition such as the Armijo rule !! $$ f(x_k + \alpha_k p_k) \le f(x_k) + c_1 \alpha_k \nabla f(x_k)^T p_k $$. class(line_search), private, allocatable :: m_lineSearch !! The line search module. logical, private :: m_useLineSearch = .true. !! Set to true if a line search should be used regardless of the !! status of m_lineSearch contains procedure, public :: get_line_search => lss_get_line_search procedure, public :: set_line_search => lss_set_line_search procedure, public :: set_default_line_search => lss_set_default procedure, public :: is_line_search_defined => & lss_is_line_search_defined procedure, public :: get_use_line_search => lss_get_use_search procedure, public :: set_use_line_search => lss_set_use_search end type type, extends(line_search_solver) :: quasi_newton_solver !! Defines a quasi-Newton solver based upon Broyden's method. The !! algorithm maintains an approximate Jacobian \(B_k\) and updates the !! iterate from the nonlinear system $F(x)=0$ using !! $$ B_k s_k = -F(x_k), \qquad x_{k+1} = x_k + s_k $$ !! with the rank-one Jacobian correction !! $$ B_{k+1} = B_k + \frac{(y_k - B_k s_k)s_k^T}{s_k^T s_k} $$ !! where \(y_k = F(x_{k+1}) - F(x_k)\). integer(int32), private :: m_jDelta = 5 !! The number of iterations that may pass between Jacobian !! calculation. contains procedure, public :: solve => qns_solve procedure, public :: get_jacobian_interval => qns_get_jac_interval procedure, public :: set_jacobian_interval => qns_set_jac_interval end type type, extends(line_search_solver) :: newton_solver !! Defines a Newton solver for systems of nonlinear equations. At each !! iteration the correction is obtained from the linearized system !! $$ J(x_k) \Delta x_k = -F(x_k), \qquad x_{k+1} = x_k + \Delta x_k $$ !! where \(J(x_k)\) is the Jacobian matrix of \(F\). contains procedure, public :: solve => ns_solve end type type, extends(equation_solver_1var) :: brent_solver !! Defines a derivative-free solver for a scalar equation \(f(x)=0\) !! based on Brent's method. The method combines bisection with secant !! and inverse quadratic interpolation to maintain a bracket and !! converge to a root. contains procedure, public :: solve => brent_solve end type type, extends(equation_solver_1var) :: newton_1var_solver !! Defines a safeguarded Newton solver for a scalar equation \(f(x)=0\). !! The iteration uses the Newton update !! $$ x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)} $$ !! and retains a bracketing interval so that the step remains bounded !! by the root bracket. contains procedure, public :: solve => newt1var_solve end type contains ! ****************************************************************************** ! LINE_SEARCH_SOLVER ! ------------------------------------------------------------------------------ subroutine lss_get_line_search(this, ls) !! Gets the line search module. class(line_search_solver), intent(in) :: this !! The [[line_search_solver]] object. class(line_search), intent(out), allocatable :: ls !! The [[line_search]] object. if (allocated(this%m_lineSearch)) & allocate(ls, source = this%m_lineSearch) end subroutine ! ---------------------- subroutine lss_set_line_search(this, ls) !! Sets the line search module. class(line_search_solver), intent(inout) :: this !! The [[line_search_solver]] object. class(line_search), intent(in) :: ls !! The [[line_search]] object. if (allocated(this%m_lineSearch)) deallocate(this%m_lineSearch) allocate(this%m_lineSearch, source = ls) end subroutine ! ------------------------------------------------------------------------------ subroutine lss_set_default(this) !! Establishes a default line_search object for the line search !! module. class(line_search_solver), intent(inout) :: this !! The [[line_search_solver]] object. type(line_search) :: ls call this%set_line_search(ls) end subroutine ! ------------------------------------------------------------------------------ pure function lss_is_line_search_defined(this) result(x) !! Tests to see if a line search module is defined. class(line_search_solver), intent(in) :: this !! The [[line_search_solver]] object. logical :: x !! Returns true if a module is defined; else, false. x = allocated(this%m_lineSearch) end function ! ------------------------------------------------------------------------------ pure function lss_get_use_search(this) result(x) !! Gets a value determining if a line-search should be employed. class(line_search_solver), intent(in) :: this !! The [[line_search_solver]] object. logical :: x !! Returns true if a line search should be used; else, false. x = this%m_useLineSearch end function ! -------------------- subroutine lss_set_use_search(this, x) !! Sets a value determining if a line-search should be employed. class(line_search_solver), intent(inout) :: this !! The [[line_search_solver]] object. logical, intent(in) :: x !! Set to true if a line search should be used; else, false. this%m_useLineSearch = x end subroutine ! ****************************************************************************** ! QUASI_NEWTON_SOLVER ! ------------------------------------------------------------------------------ subroutine qns_solve(this, fcn, x, fvec, ib, args) !! Applies the quasi-Newton's method developed by Broyden in !! conjunction with a backtracking type line search to solve N equations !! of N unknowns. !! !! See Also: !! !! - <a href="http://www.ams.org/journals/mcom/1965-19-092/S0025-5718-1965-0198670-6/S0025-5718-1965-0198670-6.pdf" !! target="_blank">Broyden's Paper</a> !! !! - <a href="https://en.wikipedia.org/wiki/Broyden%27s_method" !! target="_blank">Wikipedia</a> !! !! - <a href="http://numerical.recipes/ target="_blank">Numerical Recipes</a> class(quasi_newton_solver), intent(inout) :: this !! The [[quasi_newton_solver]] object. class(vecfcn_helper), intent(in) :: fcn !! The [[vecfcn_helper]] object containing the equations to solve. real(real64), intent(inout), dimension(:) :: x !! On input, an N-element array containing an initial estimate to !! the solution. On output, the updated solution estimate. N is !! the number of variables. real(real64), intent(out), dimension(:) :: fvec !! An N-element array that, on output, will contain the values of !! each equation as evaluated at the variable values given in x. type(iteration_behavior), optional :: ib !! An optional output, that if provided, allows the caller to !! obtain iteration performance statistics. class(*), intent(inout), optional :: args !! An optional argument to allow the user to communicate with !! fcn. ! Parameters real(real64), parameter :: zero = 0.0d0 real(real64), parameter :: half = 0.5d0 real(real64), parameter :: one = 1.0d0 real(real64), parameter :: factor = 1.0d2 ! Local Variables logical :: restart, xcnvrg, fcnvrg, gcnvrg, check integer(int32) :: i, neqn, nvar, flag, neval, iter, & maxeval, jcount, njac real(real64), allocatable, dimension(:) :: dx, df, fvold, & xold, s real(real64), allocatable, dimension(:,:) :: q, r, b real(real64) :: test, f, fold, temp, ftol, xtol, gtol, & stpmax, x2, xnorm, fnorm type(iteration_behavior) :: lib class(line_search), allocatable :: ls ! Initialization restart = .true. xcnvrg = .false. fcnvrg = .false. gcnvrg = .false. neqn = fcn%get_equation_count() nvar = fcn%get_variable_count() neval = 0 iter = 0 njac = 0 ftol = this%get_fcn_tolerance() xtol = this%get_var_tolerance() gtol = this%get_gradient_tolerance() maxeval = this%get_max_fcn_evals() if (present(ib)) then ib%iter_count = iter ib%fcn_count = neval ib%jacobian_count = njac ib%gradient_count = 0 ib%converge_on_fcn = fcnvrg ib%converge_on_chng = xcnvrg ib%converge_on_zero_diff = gcnvrg end if if (this%get_use_line_search()) then if (.not.this%is_line_search_defined()) & call this%set_default_line_search() call this%get_line_search(ls) end if ! Input Check if (.not.fcn%is_fcn_defined()) error stop NL_UNDEFINED_FUNCTION_ERROR if (nvar /= neqn) error stop NL_INVALID_INPUT_ERROR flag = 0 if (size(x) /= nvar) then flag = 3 else if (size(fvec) /= neqn) then flag = 4 end if if (flag /= 0) error stop flag ! Local Memory Allocation allocate( & b(neqn, nvar), & df(neqn), & fvold(neqn), & xold(nvar), & dx(nvar), & s(neqn) & ) ! Test to see if the initial guess is a root call fcn%fcn(x, fvec, args) f = half * dot_product(fvec, fvec) neval = neval + 1 test = zero do i = 1, neqn test = max(abs(fvec(i)), test) end do if (test < ftol) then fcnvrg = .true. end if ! Process flag = 0 ! Used to check for convergence errors if (.not.fcnvrg) then ! Determine the maximum line search step stpmax = factor * max(norm2(x), real(nvar, real64)) ! Main Iteration Loop do ! Update the iteration counter iter = iter + 1 ! Compute or update the Jacobian if (restart) then ! Compute the Jacobian call fcn%jacobian(x, b, fvec, args = args) njac = njac + 1 ! Compute the QR factorization, and form Q & R call qr_factor(b, q = q, r = r) ! Reset the Jacobian iteration counter jcount = 0 else ! Apply the rank 1 update to Q and R df = fvec - fvold dx = x - xold x2 = dot_product(dx, dx) ! Compute S = ALPHA * (DF - B * DX) s = (df - matmul(b, dx)) call recip_mult_array(x2, s) ! Compute the new Q and R matrices for the rank1 update: ! B' = B + ALPHA * S * DX**T call rank1_update(one, s, dx, b) call qr_rank1_update(q, r, s, dx) ! Increment the counter tracking how many iterations have ! passed since the last Jacobian recalculation jcount = jcount + 1 end if ! Compute GRAD = B**T * F, store in DX call mtx_mult(.true., one, b, fvec, zero, dx) ! Store FVEC and X xold = x fvold = fvec fold = f ! Solve the linear system: B * DX = -F for DX noting that ! B = Q * R. As such, form -Q**T * F, and store in DF call mtx_mult(.true., -one, q, fvec, zero, df) ! Now we have R * DX = -Q**T * F, and since R is upper ! triangular, the solution is readily computed. The solution ! will be stored in the first NVAR elements of DF df(1:nvar) = solve_triangular_system(.true., .false., .true., & r, df(1:nvar)) ! Ensure the new solution estimate is heading in a sensible ! direction. If not, it is likely time to update the Jacobian temp = dot_product(dx, df(1:nvar)) if (temp >= zero) then restart = .true. if (this%get_print_status()) then call print_status(iter, neval, njac, xnorm, fnorm) end if cycle end if ! Apply the line search if needed if (this%get_use_line_search()) then ! Define the step length for the line search temp = dot_product(df(1:nvar), df(1:nvar)) if (temp > stpmax) df(1:nvar) = df(1:nvar) * (stpmax / temp) ! Apply the line search call limit_search_vector(df(1:nvar), stpmax) call ls%search(fcn, xold, dx, df(1:nvar), x, fvec, fold, & f, lib, args = args) neval = neval + lib%fcn_count else ! No line search - just update the solution estimate x = x + df(1:nvar) call fcn%fcn(x, fvec, args) f = half * dot_product(fvec, fvec) neval = neval + 1 end if ! Test for convergence if (lib%converge_on_zero_diff .and. & this%get_use_line_search()) then call test_convergence(x, xold, fvec, dx, .true., xtol, & ftol, gtol, check, xcnvrg, fcnvrg, gcnvrg, xnorm, fnorm) else call test_convergence(x, xold, fvec, dx, .false., xtol, & ftol, gtol, check, xcnvrg, fcnvrg, gcnvrg, xnorm, fnorm) end if if (.not.check) then ! The solution did not converge, figure out why if (gcnvrg) then ! The slope of the gradient is sufficiently close to ! zero to cause issue. if (restart) then ! We've already tried recalculating a new Jacobian, ! issue a warning error stop NL_SPURIOUS_CONVERGENCE_ERROR else ! Try computing a new Jacobian restart = .true. end if else ! We have not converged, but we're not stuck with a ! zero slope gradient vector either. Go ahead and ! continue the iteration process without recomputing ! the Jacobian - unless the user dictates a ! recaclulation. if (jcount >= this%m_jDelta) then restart = .true. else restart = .false. end if end if else ! The solution has converged. It's OK to exit exit end if ! Print status if (this%get_print_status()) then call print_status(iter, neval, njac, xnorm, fnorm) end if ! Ensure we haven't made too many function evaluations if (neval >= maxeval) then flag = 1 exit end if end do end if ! Report out iteration statistics if (present(ib)) then ib%iter_count = iter ib%fcn_count = neval ib%jacobian_count = njac ib%gradient_count = 0 ib%converge_on_fcn = fcnvrg ib%converge_on_chng = xcnvrg ib%converge_on_zero_diff = gcnvrg end if ! Check for convergence issues if (flag /= 0) then error stop NL_CONVERGENCE_ERROR end if end subroutine ! ------------------------------------------------------------------------------ pure function qns_get_jac_interval(this) result(n) !! Gets the number of iterations that may pass before forcing a !! recalculation of the Jacobian matrix. class(quasi_newton_solver), intent(in) :: this !! The [[quasi_newton_solver]] object. integer(int32) :: n !! The number of iterations. n = this%m_jDelta end function ! -------------------- subroutine qns_set_jac_interval(this, n) !! Sets the number of iterations that may pass before forcing a !! recalculation of the Jacobian matrix. class(quasi_newton_solver), intent(inout) :: this !! The [[quasi_newton_solver]] object. integer(int32), intent(in) :: n !! The number of iterations. this%m_jDelta = n end subroutine ! ****************************************************************************** ! NEWTON_SOLVER ! ------------------------------------------------------------------------------ subroutine ns_solve(this, fcn, x, fvec, ib, args) !! Applies Newton's method in conjunction with a backtracking type !! line search to solve N equations of N unknowns. class(newton_solver), intent(inout) :: this !! The [[newton_solver]] object. class(vecfcn_helper), intent(in) :: fcn !! The [[vecfcn_helper]] object containing the equations to solve. real(real64), intent(inout), dimension(:) :: x !! On input, an N-element array containing an initial estimate to !! the solution. On output, the updated solution estimate. N is !! the number of variables. real(real64), intent(out), dimension(:) :: fvec !! An N-element array that, on output, will contain the values of !! each equation as evaluated at the variable values given in x. type(iteration_behavior), optional :: ib !! An optional output, that if provided, allows the caller to !! obtain iteration performance statistics. class(*), intent(inout), optional :: args !! An optional argument to allow the user to communicate with fcn. ! Parameters real(real64), parameter :: zero = 0.0d0 real(real64), parameter :: half = 0.5d0 real(real64), parameter :: one = 1.0d0 real(real64), parameter :: mintol = 1.0d-12 real(real64), parameter :: factor = 1.0d2 ! Local Variables logical :: check, xcnvrg, fcnvrg, gcnvrg integer(int32) :: i, neqn, nvar, flag, neval, iter, maxeval, njac integer(int32), allocatable, dimension(:) :: ipvt real(real64), allocatable, dimension(:) :: dir, grad, xold real(real64), allocatable, dimension(:,:) :: jac, lu real(real64) :: ftol, xtol, gtol, f, fold, stpmax, xnorm, fnorm, temp, test type(iteration_behavior) :: lib class(line_search), allocatable :: ls ! Initialization xcnvrg = .false. fcnvrg = .false. gcnvrg = .false. neqn = fcn%get_equation_count() nvar = fcn%get_variable_count() neval = 0 iter = 0 njac = 0 ftol = this%get_fcn_tolerance() xtol = this%get_var_tolerance() gtol = this%get_gradient_tolerance() maxeval = this%get_max_fcn_evals() if (present(ib)) then ib%iter_count = iter ib%fcn_count = neval ib%jacobian_count = njac ib%gradient_count = 0 ib%converge_on_fcn = fcnvrg ib%converge_on_chng = xcnvrg ib%converge_on_zero_diff = gcnvrg end if if (this%get_use_line_search()) then if (.not.this%is_line_search_defined()) & call this%set_default_line_search() call this%get_line_search(ls) end if ! Input Checking if (.not.fcn%is_fcn_defined()) error stop NL_UNDEFINED_FUNCTION_ERROR if (nvar /= neqn) error stop NL_INVALID_INPUT_ERROR flag = 0 if (size(x) /= nvar) then flag = 3 else if (size(fvec) /= neqn) then flag = 4 end if if (flag /= 0) error stop flag ! Local Memory Allocation allocate( & dir(nvar), & grad(nvar), & xold(nvar), & jac(nvar, neqn) & ) call fcn%jacobian(x, jac, fv = fvec, args = args) ! Test to see if the initial guess is a root call fcn%fcn(x, fvec, args) f = half * dot_product(fvec, fvec) neval = neval + 1 test = zero do i = 1, neqn test = max(abs(fvec(i)), test) end do if (test < ftol) then fcnvrg = .true. end if ! Process flag = 0 ! Used to check for convergence errors if (.not.fcnvrg) then ! Compute the maximum step size for the line search process stpmax = factor * max(norm2(x), real(nvar, real64)) ! Main Iteration Loop do ! Increment the iteration counter iter = iter + 1 ! Compute the Jacobian call fcn%jacobian(x, jac, fvec, args = args) njac = njac + 1 ! Compute the gradient do i = 1, nvar grad(i) = dot_product(jac(:,i), fvec) end do ! Compute the LU factorization of the Jacobian call lu_factor(jac, ipvt = ipvt, lu = lu) ! Store previous iteration values xold = x fold = f ! Solve the linear system of equations dir = solve_lu(lu, ipvt, -fvec) ! Apply the line search if needed if (this%get_use_line_search()) then ! Define the step length for the line search temp = dot_product(dir, dir) if (temp > stpmax) dir = dir * (stpmax / temp) ! Apply the line search call limit_search_vector(dir, stpmax) call ls%search(fcn, xold, grad, dir, x, fvec, & fold, f, lib, args = args) neval = neval + lib%fcn_count else ! No line search - just update the solution estimate x = x + dir call fcn%fcn(x, fvec, args) f = half * dot_product(fvec, fvec) neval = neval + 1 end if ! Check for convergence call test_convergence(x, xold, fvec, grad, .true., xtol, & ftol, gtol, check, xcnvrg, fcnvrg, gcnvrg, xnorm, fnorm) if (check) then ! The solution has converged exit else if (gcnvrg) then ! The solution appears to have settled at a point where ! the gradient has a zero slope error stop NL_SPURIOUS_CONVERGENCE_ERROR end if ! Print status if (this%get_print_status()) then call print_status(iter, neval, njac, xnorm, fnorm) end if ! Ensure we haven't made too many function evaluations if (neval >= maxeval) then flag = 1 exit end if end do end if ! Report out iteration statistics if (present(ib)) then ib%iter_count = iter ib%fcn_count = neval ib%jacobian_count = njac ib%gradient_count = 0 ib%converge_on_fcn = fcnvrg ib%converge_on_chng = xcnvrg ib%converge_on_zero_diff = gcnvrg end if ! Check for convergence issues if (flag /= 0) then error stop NL_CONVERGENCE_ERROR end if end subroutine ! ****************************************************************************** ! BRENT_SOLVER ! ------------------------------------------------------------------------------ subroutine brent_solve(this, fcn, x, lim, f, ib, args) !! Solves an equation of one variable using Brent's method. !! !! See Also !! !! - <a href="https://en.wikipedia.org/wiki/Brent%27s_method" !! target="_blank">Wikipedia</a> !! !! - <a href="http://numerical.recipes/ target="_blank">Numerical Recipes</a> !! !! - <a href="https://maths-people.anu.edu.au/~brent/pub/pub011.html" !! target="_blank">R.P. Brent, "Algorithms for Minimization without !! Derivatives."</a> class(brent_solver), intent(inout) :: this !! The [[brent_solver]] object. class(fcn1var_helper), intent(in) :: fcn !! The [[fcn1var_helper]] object containing the equation to solve. real(real64), intent(inout) :: x !! A parameter used to return the solution. Notice, any input !! value will be ignored as this routine relies upon the search !! limits in lim to provide a starting point. type(value_pair), intent(in) :: lim !! A [[value_pair]] object defining the search limits. real(real64), intent(out), optional :: f !! An optional parameter used to return the function residual as !! computed at x. type(iteration_behavior), optional :: ib !! An optional output, that if provided, allows the caller to !! obtain iteration performance statistics. class(*), intent(inout), optional :: args !! An optional argument to allow the user to communicate with fcn. ! Parameters real(real64), parameter :: zero = 0.0d0 real(real64), parameter :: half = 0.5d0 real(real64), parameter :: one = 1.0d0 real(real64), parameter :: two = 2.0d0 real(real64), parameter :: three = 3.0d0 ! Local Variables logical :: fcnvrg, xcnvrg integer(int32) :: neval, maxeval, flag, iter real(real64) :: ftol, xtol, a, b, c, fa, fb, fc, p, q, r, s, xm, e, d, & mn1, mn2, eps, tol1, temp ! Initialization fcnvrg = .false. xcnvrg = .false. x = zero a = min(lim%x1, lim%x2) b = max(lim%x1, lim%x2) neval = 0 iter = 0 eps = epsilon(eps) ftol = this%get_fcn_tolerance() xtol = this%get_var_tolerance() maxeval = this%get_max_fcn_evals() if (present(f)) f = zero if (present(ib)) then ib%iter_count = iter ib%fcn_count = neval ib%jacobian_count = 0 ib%gradient_count = 0 ib%converge_on_fcn = fcnvrg ib%converge_on_chng = xcnvrg ib%converge_on_zero_diff = .false. end if ! Input Check if (.not.fcn%is_fcn_defined()) error stop NL_UNDEFINED_FUNCTION_ERROR if (abs(a - b) < eps) error stop NL_INVALID_INPUT_ERROR ! Process flag = 0 fa = fcn%fcn(a, args) fb = fcn%fcn(b, args) neval = 2 fc = fb do ! Increment the iteration counter iter = iter + 1 ! Adjust the bounding interval if ((fb > zero .and. fc >= zero) .or. & (fb < zero .and. fc < zero)) then c = a fc = fa d = b - a e = d end if if (abs(fc) < abs(fb)) then a = b b = c c = a fa = fb fb = fc fc = fa end if ! Convergence Check tol1 = two * eps * abs(b) + half * xtol xm = half * (c - b) if (abs(fb) < ftol) then x = b fcnvrg = .true. exit end if if (abs(xm) <= tol1) then x = b xcnvrg = .true. exit end if ! Actual Method if (abs(e) >= tol1 .and. abs(fa) > abs(fb)) then ! Attempt the inverse quadratic interpolation to determine ! the root s = fb / fa if (abs(a - c) < eps) then ! a == c p = two * xm * s q = one - s else q = fa / fc r = fb / fc p = s * (two * xm * q * (q - r) - (b - a) * (r - one)) q = (q - one) * (r - one) * (s - one) end if ! Ensure we're within bounds if (p > zero) q = -q p = abs(p) mn1 = three * xm * q - abs(tol1 * q) mn2 = abs(e * q) if (mn1 < mn2) then temp = mn1 else temp = mn2 end if if (two * p < temp) then ! Accept the interpolation e = d d = p / q else ! The interpolation failed, use bisection d = xm e = d end if else ! The bounds are decreasing too slowly, use bisection d = xm e = d end if ! Move the last best guess to the lower limit parameter (A) a = b fa = fb if (abs(d) > tol1) then b = b + d else b = b + sign(tol1, xm) end if fb = fcn%fcn(b, args) neval = neval + 1 ! Print iteration status if (this%get_print_status()) then call print_status(iter, neval, 0, xm, fb) end if ! Ensure we haven't made too many function evaluations if (neval >= maxeval) then flag = 1 exit end if end do ! Report out iteration statistics and other optional outputs if (present(f)) f = fb if (present(ib)) then ib%iter_count = iter ib%fcn_count = neval ib%jacobian_count = 0 ib%gradient_count = 0 ib%converge_on_fcn = fcnvrg ib%converge_on_chng = xcnvrg ib%converge_on_zero_diff = .false. end if ! Check for convergence issues if (flag /= 0) then error stop NL_CONVERGENCE_ERROR end if end subroutine ! ****************************************************************************** ! NEWTON_1VAR_SOLVER ! ------------------------------------------------------------------------------ subroutine newt1var_solve(this, fcn, x, lim, f, ib, args) !! Solves an equation of one variable using Newton's method. class(newton_1var_solver), intent(inout) :: this !! The [[newton_1var_solver]] object. class(fcn1var_helper), intent(in) :: fcn !! The [[fcn1var_helper]] object containing the equation to solve. real(real64), intent(inout) :: x !! A parameter used to return the solution. Notice, any input !! value will be ignored as this routine relies upon the search !! limits in lim to provide a starting point. type(value_pair), intent(in) :: lim !! A value_pair object defining the search limits. real(real64), intent(out), optional :: f !! An optional parameter used to return the function residual as !! computed at x. type(iteration_behavior), optional :: ib !! An optional output, that if provided, allows the caller to !! obtain iteration performance statistics. class(*), intent(inout), optional :: args !! An optional argument to allow the user to communicate with fcn. ! Parameters real(real64), parameter :: zero = 0.0d0 real(real64), parameter :: p5 = 0.5d0 real(real64), parameter :: two = 2.0d0 ! Local Variables logical :: fcnvrg, xcnvrg, dcnvrg integer(int32) :: neval, ndiff, maxeval, flag, iter real(real64) :: ftol, xtol, dtol, xh, xl, fh, fl, x1, x2, eps, dxold, & dx, df, temp, ff ! Initialization fcnvrg = .false. xcnvrg = .false. dcnvrg = .false. neval = 0 ndiff = 0 iter = 0 ftol = this%get_fcn_tolerance() xtol = this%get_var_tolerance() dtol = this%get_diff_tolerance() maxeval = this%get_max_fcn_evals() if (present(f)) f = zero if (present(ib)) then ib%iter_count = iter ib%fcn_count = neval ib%jacobian_count = ndiff ib%gradient_count = 0 ib%converge_on_fcn = fcnvrg ib%converge_on_chng = xcnvrg ib%converge_on_zero_diff = dcnvrg end if x1 = min(lim%x1, lim%x2) x2 = max(lim%x1, lim%x2) eps = epsilon(eps) ! Input Check if (.not.fcn%is_fcn_defined()) error stop NL_UNDEFINED_FUNCTION_ERROR if (abs(x1 - x2) < eps) error stop NL_INVALID_INPUT_ERROR ! See if the root is one of the end points flag = 0 fl = fcn%fcn(x1, args) fh = fcn%fcn(x2, args) neval = 2 if (abs(fl) < ftol) then x = x1 if (present(f)) f = fl if (present(ib)) then ib%converge_on_fcn = .true. ib%fcn_count = 2 end if return end if if (abs(fh) < ftol) then x = x2 if (present(f)) f = fh if (present(ib)) then ib%converge_on_fcn = .true. ib%fcn_count = 2 end if return end if ! Process if (fl < zero) then xl = x1 xh = x2 else xl = x2 xh = x1 end if x = p5 * (x1 + x2) dxold = abs(x2 - x1) dx = dxold ff = fcn%fcn(x, args) df = fcn%diff(x, f = ff, args = args) neval = neval + 1 ndiff = ndiff + 1 do ! Increment the iteration counter iter = iter + 1 ! Bisect if the Newton step went out of range, or if the rate ! of change was too slow if ((((x - xh) * df - ff) * ((x - xl) * df - ff) > zero) .or. & (abs(two * ff) > abs(dxold * df))) & then ! Bisection dxold = dx dx = p5 * (xh - xl) x = xl + dx if (abs(xl - x) < xtol) then ! Convergence as the change in root is within tolerance xcnvrg = .true. exit end if else ! Newton's Method dxold = dx dx = ff / df temp = x x = x - dx if (abs(temp - x) < xtol) then ! Convergence as the change in root is within tolerance xcnvrg = .true. exit end if end if ! Update function values ff = fcn%fcn(x, args) df = fcn%diff(x, f = ff, args = args) neval = neval + 1 ndiff = ndiff + 1 ! Check for convergence if (abs(ff) < ftol) then fcnvrg = .true. exit end if if (abs(dx) < xtol) then xcnvrg = .true. exit end if if (abs(df) < dtol) then dcnvrg = .true. exit end if ! Update the bracket on the root if (ff < zero) then xl = x else xh = x end if ! Print status if (this%get_print_status()) then call print_status(iter, neval, ndiff, dx, ff) end if ! Ensure we haven't made too many function evaluations if (neval >= maxeval) then flag = 1 exit end if end do ! Ensure the function value is current with the estimate of the root if (present(f)) then f = fcn%fcn(x, args) neval = neval + 1 end if ! Report out iteration statistics and other optional outputs if (present(f)) f = ff if (present(ib)) then ib%iter_count = iter ib%fcn_count = neval ib%jacobian_count = ndiff ib%gradient_count = 0 ib%converge_on_fcn = fcnvrg ib%converge_on_chng = xcnvrg ib%converge_on_zero_diff = dcnvrg end if ! Check for convergence issues if (flag /= 0) then error stop NL_CONVERGENCE_ERROR end if end subroutine end module